
% Function that computes sample moments

function mhat = mhat(ln_y1,ln_y2,ln_a,ln_ri,ln_rc)

    % Matrix to extract covariance from variance-covariance matrix
    M           = [0,0;1,0];

    % Sample moments
    mu_y1       = mean(ln_y1);
    mu_y2       = mean(ln_y2);
    mu_a        = mean(ln_a);
    mu_ri       = mean(ln_ri);
    mu_rc       = mean(ln_rc);
    var_y1      = var(ln_y1);
    var_y2      = var(ln_y2);
    var_a       = var(ln_a);
    var_ri      = var(ln_ri);
    var_rc      = var(ln_rc);
    cov_y1y2    = sum(cov(ln_y1,ln_y2).*M,'all');
    cov_y1a     = sum(cov(ln_y1,ln_a).*M,'all');
    cov_y1ri    = sum(cov(ln_y1,ln_ri).*M,'all');
    cov_y1rc    = sum(cov(ln_y1,ln_rc).*M,'all');
    cov_y2a     = sum(cov(ln_y2,ln_a).*M,'all');
    cov_y2ri    = sum(cov(ln_y2,ln_ri).*M,'all');
    cov_y2rc    = sum(cov(ln_y2,ln_rc).*M,'all');
    cov_ari     = sum(cov(ln_a,ln_ri).*M,'all');
    cov_arc     = sum(cov(ln_a,ln_rc).*M,'all');
    cov_rirc    = sum(cov(ln_ri,ln_rc).*M,'all');

    % Vector of sample moments
    mhat        = [mu_y1;mu_y2;mu_a;mu_ri;mu_rc;var_y1;cov_y1y2;...
                cov_y1a;cov_y1ri;cov_y1rc;var_y2;cov_y2a;cov_y2ri;...
                cov_y2rc;var_a;cov_ari;cov_arc;var_ri;cov_rirc;var_rc];
            